MODULE md1_main

   USE md1_global_consts
   USE md1_global_types
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: md1_initial,                        &
             md1_final,                          &
             md1_write_multipole_integrals,      &
             md1_get_AO_FE_matrix_elements,      &
             md1_get_FE_charge_overlap

   TYPE(md1_basis), SAVE, POINTER :: basis

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE md1_initial(filename,nbas)

      USE md1_basis_set, ONLY: md1_init_basis_set, md1_get_basis_set

      IMPLICIT NONE
      CHARACTER(*),  INTENT(IN) :: filename    
      INTEGER(INTK), INTENT(IN) :: nbas

      WRITE(IO1,*)
      WRITE(IO1,*), "-----------------------------------------"
      WRITE(IO1,*), "|         MD1 initialization            |"
      WRITE(IO1,*), "-----------------------------------------"
      WRITE(IO1,*)
      NULLIFY(basis)
      CALL md1_init_basis_set(filename)
      CALL md1_get_basis_set(basis)
      IF (nbas /= basis%nbas) STOP 'md1_initial: nbas error!'

   END SUBROUTINE md1_initial

!-------------------------------------------------------------------------------

   SUBROUTINE md1_final

      USE md1_basis_set,   ONLY: md1_free_basis_set
      USE md1_shell_pairs, ONLY: md1_free_shell_pairs

      IMPLICIT NONE
      CALL md1_free_basis_set
      CALL md1_free_shell_pairs
      NULLIFY(basis)

   END SUBROUTINE md1_final

!-------------------------------------------------------------------------------

   SUBROUTINE md1_write_multipole_integrals(LMAX)

      USE md1_shell_pairs,    ONLY: md1_get_shell_pairs
      USE md1_multipole_ints, ONLY: md1_init_multipole_ints,            &
                                    md1_free_multipole_ints,            &
                                    md1_build_and_write_sph_mpoles

      IMPLICIT NONE  
      INTEGER(INTK), INTENT(IN) :: LMAX

      TYPE(md1_sh_pairs), POINTER :: sh_pairs(:)
      INTEGER(INTK) :: nmoms

      ! Get list of non-vanishing shell pairs
      CALL md1_get_shell_pairs(basis,sh_pairs)

      ! Build and store GTO multipole moments on disk
      WRITE(IO1,*), "MD1 multipole integral LMAX =", LMAX

      ! Open multipole file for interface with fmmlib
      OPEN(UNIT=IO2, FILE='multipoles.md1', STATUS='REPLACE',  &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IO2)                   
                                    
      ! Write overlap matrix to disk for diagnostics
      OPEN(UNIT=IO3, FILE="overlap_matrix.md1", STATUS='REPLACE',  &
           ACCESS='SEQUENTIAL', FORM='FORMATTED')
      REWIND(IO3)

      CALL md1_init_multipole_ints(basis,LMAX)
      CALL md1_build_and_write_sph_mpoles(basis,LMAX,sh_pairs,nmoms)
      CALL md1_free_multipole_ints

      ! Mark end of file with negative L-index
      WRITE(IO2) 0,-1,0,0,0, 0d0,0d0,0d0, 0d0
      CLOSE(UNIT=IO2, STATUS='KEEP')
      WRITE(IO3,'(3I8,F25.15)') 0,-1,0, 0d0
      CLOSE(UNIT=IO3, STATUS='KEEP')

      ! Write header file with basic constants and dimensions
      OPEN(UNIT=IO2, FILE='multipoles_header.md1', STATUS='REPLACE',   &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      WRITE(IO2) LMAX, basis%nbas, nmoms
      CLOSE(UNIT=IO2, STATUS='KEEP')

      NULLIFY(sh_pairs)

   END SUBROUTINE md1_write_multipole_integrals

!-------------------------------------------------------------------------------

   SUBROUTINE md1_get_AO_FE_matrix_elements(grid,Q_matrix)

      USE md1_shell_pairs,   ONLY: md1_get_shell_pairs
      USE md1_FE_AO_overlap, ONLY: md1_get_FE_AO_overlaps

      IMPLICIT NONE  
      TYPE(md1_FE_grid), INTENT(IN)  :: grid
      REAL(REALK),       INTENT(OUT) :: Q_matrix(:,:)

      TYPE(md1_sh_pairs), POINTER :: sh_pairs(:)
      INTEGER(INTK) :: i

      ! Get list of non-vanishing shell pairs
      CALL md1_get_shell_pairs(basis,sh_pairs)
      ! Overlap FE field with all AO pairs
      CALL md1_get_FE_AO_overlaps(basis,sh_pairs,grid,Q_matrix)

      NULLIFY(sh_pairs)

   END SUBROUTINE md1_get_AO_FE_matrix_elements

!-------------------------------------------------------------------------------

   SUBROUTINE md1_get_FE_charge_overlap(n,charges,centres,grid,energies)

      USE md1_LIP_coeffs,    ONLY: md1_init_LIP_coeffs, md1_free_LIP_coeffs
      USE md1_FE_AO_overlap, ONLY: md1_overlap_AO_FE

      IMPLICIT NONE
      INTEGER(INTK),     INTENT(IN)  :: n
      REAL(REALK),       INTENT(IN)  :: charges(n)
      REAL(REALK),       INTENT(IN)  :: centres(3,n)
      TYPE(md1_FE_grid), INTENT(IN)  :: grid
      REAL(REALK),       INTENT(OUT) :: energies(n,n)

      REAL(REALK), PARAMETER :: nuc_exponent = 1d6
      REAL(REALK), PARAMETER :: normalization = (SQRT(PI/1d6)**3)
      TYPE(md1_sh_pairs)   :: pseudo_sh
      TYPE(md1_prim_batch) :: batch(1)
      INTEGER(INTK) :: k

      CALL md1_init_LIP_coeffs(grid%order,grid%h)

      ! Map charges onto AO pairs data structures to re-use old code
      DO k = 1, SIZE(charges)
         ! "AO shell pair"
         pseudo_sh%P%lo = k
         pseudo_sh%P%hi = k
         pseudo_sh%Q%lo = k
         pseudo_sh%Q%hi = k
         pseudo_sh%xyz = centres(:,k)
         pseudo_sh%I = k
         pseudo_sh%J = k
         pseudo_sh%Iang = 0
         pseudo_sh%Jang = 0
         ! "primitive AO pairs"
         batch(1)%P = centres(:,k)
         batch(1)%PA = centres(:,k)
         batch(1)%PB = centres(:,k)
         ! Model point charge as a very compact Gaussian
         batch(1)%exp = nuc_exponent
         batch(1)%CAB = charges(k) / normalization
         batch(1)%KAB = 1d0
         batch(1)%fmm_extent = zero
         ! Call code for AO/FE overlap
         CALL md1_overlap_AO_FE(basis%LM,pseudo_sh,batch,grid,energies)
      END DO

      CALL md1_free_LIP_coeffs

   END SUBROUTINE md1_get_FE_charge_overlap

!-------------------------------------------------------------------------------

END MODULE md1_main
